Draft version February 22, 2012 

Preprint typeset using 1^1^^ style emulatcapj v. 02/07/07 



T— I 
O 
(N 

X) 
D 

u: 
^' 

o 

(N 
> 

O 



A DENSITY INDEPENDENT FORMULATION OF SMOOTHED PARTICLE HYDRODYNAMICS 

Takayuki R. Saitoh^ & Junichiro Making^ 

Draft version February 22, 2012 

ABSTRACT 

In the standard formulation of the smoothed particle hydrodynamics (SPH), it is assumed that the 
local density distribution is differentiable. This assumption is used to derive the spatial derivatives of 
other quantities. However, this assumption breaks down at the contact discontinuity, which appears 
often in simulations of astronomical objects. At the contact discontinuity, the density of the low- 
density side is overestimated while that of the high-density side is underestimated. As a result, the 
pressure of the low (high) density side is over (under) estimated. Thus, unphysical repulsive force 
appears at the contact discontinuity, resulting in the effective surface tension. This effective surface 
tension suppresses instabilities such as the Kelvin-Hclmholtz and Rayleigh- Taylor instabilities. In this 
paper, we present a new formulation of SPH, which does not require the differentiability of density 
and thus can handle contact discontinuity without numerical problems. The results of standard tests 
such as the shock tube, Kelvin-Hclmholtz and Rayleigh- Taylor instabilities, and the blob tests are all 
very favorable to our new formulation. We conclude that our new formulation solved practically all 
known difficulties of the standard SPH, without introducing additional numerical diffusion or breaking 
the exact force symmetry or energy conservation. 

Subject headings: galaxies: evolution — galaxies:ISM — methods:numerical 



1. INTRODUCTION 

Smoothed particle hydrodynamics (SPH) is a lagrange 
scheme to solve the evolution of f luid using p arti- 
cles. It was originally intro duced by iLucvl ()1977t ) and 
iGingold fc Monaghanl (|1977| ) and has been widely used 
in the field of the computa t ional astrophysi cs ([Monaghanl 
fl99l [20051: rRosswog|[200l lSpringe]|[2010bl) . It is becom- 
ing popu lar in hydrodyna mical simulations in engineer- 
ing re.g.. lLiu fc Lh][2n0^ 

Rccentlv. lAgertz et al.l (I2007D reported the results of 
comparison of SPH and Euler schemes (grid methods). 
Their main result is that SPH suppresses the Kclvin- 
Helmholtz inst a bility. This has been pointed out by 
lOkamoto etahl ([2003[) . The reason of this problem is 
that in the standard SPH the smoothed density is used 
to obtain other physical quantities. The estimated den- 
sity of particles near the contact discontinuity has 0(1) 
error, irrespective of the numerical resolution. This large 
error causes s i milar ly large error in the pressure (see 
lAgertz et al.l ([2007() noted that there were fundamental 
differences between SPH and grid methods. 

There have been several proposals to improve S PH so 
that i t can deal with the contact discontinuity. iPricd 
([20081) discussed the artificia l thermal conducti vity which 
was originally introduced bv lMonaghanl ([19971 ). The mo- 
tivation of the use of the artificial conductivity is that ev- 
ery physical quantity should be smooth in the standard 
SPH. The artificial conductivity eliminates the discon- 
tinuity in the thermal energy. As a result, the density 
near the contact discontinuity becomes smooth and the 
pressure also becomes smooth with this artificial con- 
ductivity. Thus, the Kclvin-Helmholtz instability takes 
place. At the first sight, this artificial conductivity looks 
similar to the artificial viscosity which is necessary to 
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capture shocks in SPH. However, there are two funda- 
mental differences. First, the artificial viscosity is used 
to generate the physical dissipation associated with the 
shock, while the artificial conductivity adds physically 
non-existent dissipation. One needs to fine-tune the con- 
ductivity coefficient to prevent unnecessary smoothing. 
This means that the conductivity must be non-linear. 
Second, if there is a jump in the chemical composition, 
thermal conductivity is not enough. However, whether 
the use of artificial chemical di f fusion is justified or not is 
an open question. iRead et al.l ([2010D suggested that the 
Kelvin-Helmholtz instability took place when the num- 
ber of neighbors wa s sufficiently large and th e momen- 
tum equatio of the iRitchie fc Thoma^ ([20011 ) was used. 
lAbell ([20 lit ) used the relative pressure instead of the 
absolute values of pressures in the equation of motion. 
This formulation improves the treatment of the Kelvin- 
Helmholtz instabihty, bu t breaks the Newton's third law. 
I Garcia- Senz et all (|2012| ) considered the use of the inte- 
gral form of the first derivative, which also improved the 
treatment of hydrodynamical instabilities. 

In this paper, we describe a new formulation of SPH 
which does not use the smoothed mass density. Instead, 
we use the smoothed internal energy density to obtain 
other quantities and their spatial derivatives. The rea- 
son why we adopt the energy density instead of the mass 
density is that it is the fundamental quantity of the hy- 
drodynamics. In our formulation, the pressure is calcu- 
lated without using the mass density. Thus, unphysi- 
cal jumps of pressure at the contact discontinuity dis- 
appcar. Our equation o f motion is similar to that of 
iRitchie fc Thomas! ([2001|) . although the basis of devia- 
tion is completely different. Results of various tests in- 
dicate that our formulation is highly advantageous. 

The structure of this paper is as follows. In fJ21 we 
analyze the problem of standard SPH at discontinuities. 
Our new formulation of SPH is described in fj3] and the 
comparison of the results of test calculations with the 
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new formulation and standard formulation arc shown in 
^ Summary and discussion are presented in ^ 

2. STANDARD SPH AND ITS DIFFICULTY AROUND 
DISCONTINUITIES 

In SPH, the fluid is expressed by discrete particles and 
physical quantities are approximated by kernel interpo- 
lation. In the standard formulation of SPH, the local 
density is first calculated, and then the rests of neces- 
sary physical quantities, such as the pressure gradient 
and the time derivative of the internal energy, are calcu- 
lated. Thus, the accuracy of the solution depends on the 
accuracy of the density estimate. In this section, we re- 
examine the derivation of the equation of motion in SPH 
to understand its problem. 

A physical quantity / at position r can be expressed 
as follows: 



fir) 



f{r')5{\r-r'\)dr'. 



(1) 



A smoothed value of / at position r, (/)(r), is given by 
the convolution of / and a kernel function W{r — r', h): 



{f)ir)= / fir')Wi\r-r'\,h)dr', 



(2) 



where h is the size of the kernel function and corresponds 
to the spatial resolution. This smoothing is the base of 
SPH. Here, the kernel function must satisfy the following 
three conditions: (1) it becomes the delta function in the 
limit of /i — )■ 0, (2) it is normalized as unity, and (3) it is a 
function with compact support. A cubic spline function 
is most widely used as the kernel function: 



W{\r-r'lh) 




< s < 1, 

1 < s < 2, (3) 

2 < s. 



where s — \r — r'l/h, D is the dimension, and the nor- 
malized factors a in one, two, and three dimensions are 
2/3, 10/77r, and I/tt, respectively. We first derive the 
equations of motion and energy with the constant ker- 
nel size, and then we generalized them to the individual 
kernel size. 

The first derivative of the smoothed / is given by 



(V/)(r) = J S/f{r')W{\ 



r — r'l, h)dr' . 



(4) 



By making use of the partial integral and the fact that 
the kernel function has compact support, Eq. 2] becomes 



(V/>(r)= / f{r')\/W{\r~r'\,h)dr'. 



(5) 



We need to discretize Eq. [2] to evaluate the physical 
quantities at positions of particles. To convert integral 
into summation, a volume element dr' is replaced by 
m j / pj , where nij and pj are the mass and density of the 
particle j. In addition, positions of particles i and j are 
expressed by and rj and /(r') is replaced by fj. Thus, 
the value of / at the position of particle i is 



(/)(r,)c.^m,^Ty(n„/i), 
Pj 



(6) 



the standard SPH. At this point, we do not know pj. By 
substituting p into /, we obtain 



Pi ~ ^mj H/(ry ,/i), 



(7) 



where pi = is the smoothed density at the posi- 

tion of particle i. Note that the right-hand side of Eq. [7] 
includes no unknown quantities. Thus, densities should 
be calculated first in the standard SPH. 
The equation of motion is 



d^r _ VP 



(8) 



where t is time and P is pressure. The SPH approxima- 
tion of Eq. IS] is given by 



, \Pi Pj I 



dt^ 



This form satisfies the Newton's third law. We used the 
following relation to obtain Eq. [HI 



VP ^fP\ P^ 
p \p J p^ 



(10) 



In order for Eq. [2] to be meaningful, p must be differen- 
tiable, since its derivative is used in Eq. 1101 

Finally, we derive the energy equation in the standard 
SPH. The energy equation is 

P_ , ^ 

— = V-t;, (11) 

dt p 

where u is the internal energy and v is the velocity. To 
obtain the SPH formulation of the energy equation, we 
need the SPH expression of V • t>. We use 

V{pv) ^Vpv + pV ■ V. (12) 

The SPH formulation of V • t; is given by 

PjV • ~ ^ m-iVj ■ WWinj ,h) -Vi-^ nij'VWin.i , h) 



^mjUy • VW^(r,j,/i), 



(13) 



where Vij = Vi — Vj. Therefore, the energy equation in 
the standard SPH is 

^^Y.'^A'^^r^Win,,h). (14) 

ai Pi 

Equations [71 [HI and [U close with the equation of state 
(EOS), 

P=(7-1)PM, (15) 

where 7 is the specific heat ratio. There is no need to 
solve the continuity equation in SPH since it is satisfied 
automatically. 

When we use the variable and individual kernel size, 
above equations should be modified slightly. First, the 
density evaluation equation is rewritten as 



Pi 



^mjW{nj,h,). 



(16) 



where = |ry | and = — r^. Hereafter, we call 
the SPH formulation with this type of discretization as 



This is the so-called gather int erpretation of the summa- 
tion (jHernguist &: Kat3ll989l ). In equations of motion 
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and energy, the gather-an d-scatter interpretation is used 
(|Hernquist fc Katdll989l ). Thus, Eqs. Eland [H become 



and 



dui 
It 



(17) 



(18) 



where VW{rij , hi) is replaced by VWij = 
0.5[yW{rij,hi) + VW{rij,hj)] so that the equa- 
tion of motion can satisfy the Newton's third law. It is 
also possible to use \7Wij = VW[rij ,0.5{hi + hj)]. We 
adopt the first form throughout this paper. 

In the derivation of the standard SPH discretization, 
the differentiability of p is used for both the equation of 
motion and the energy equation. However, p is discon- 
tinuous at the contact discontinuity. In the following, 
we illustrate the consequence of the discontinuity of the 
density. 

In Figure (U we show the values of density and pres- 
sure around a contact discontinuity evaluated by the 
standard formulation of SPH. Equation [16] is used and 
P = (7 — l)pu. To set up this contact discontinuity, we 
place particles on a regular grid in three dimensions and 
set p = 1 for a; < 0.5 and p = 0.125 for x > 0.5. We used 
equal-mass particles in the first two configurations. In 
these two setups, positions of particles in the less dense 
region is determined by taking the center of mass of the 
eight particles in the cube of the particle separation. In 
the last configuration, we adopted the equal separation 
for both regions, which means that the mass of parti- 
cles in the less dense region is 1/8 of that of particles 
in the dense region. The internal energy was set to 1.5 
{x < 0.5) and 12 {x > 0.5), and the specific heat ratio 
was 5/3. Velocities of particles were set to zero. The 
kernel size is determined to keep the neighbor number, 
A'^nb, to the range 32 ± 2, in the first and the last tests. 
In the second test a constant h fixed to twice the particle 
separation in the less dense region is used. 

The top panels show the distribution of particles. The 
panels in the second row show the SPH density. Though 
the initial setup has the discontinuity at x = 0.5, it 
is smoothed by the kernel. As a result, SPH density 
of particles next to the discontinuity has very large er- 
rors, as shown in the 3rd row. This large error in 
the density causes similarly large error in the pressure 
(4th row). The pressure of particles at the end of the 
low-density region is grossly overestimated, while that 
at the end of the high-density region is underestimated 
only modestly. This non-symmetric error in the pres- 
sure is the origin of the repulsive force at the contact 
discontin uity, as has been pointed out in previous stud- 
ies (e.g., iRitchie fc Thomasll200ll lOkamoto et al.ll2003l: 
lAgertz et al.ll2007( ). This large error in the pressure also 
exists in both of the constant kernel size case (the mid- 
dle column) and the equal-separation case (the right col- 
umn). 

Consider the following density and pressure distribu- 
tion: 

p=h (19) 

\p2 x<0, 



and 



Obviously, we have 



P = Po- 



{p){x) ^ t^' , forx->0. 



and therefore. 



lim(P)(.).^Po, 

x->-+0 Zpi 

lim(P)(.) = ^Po. 

x^-a ip2 



(20) 
(21) 

(22) 
(23) 



Thus, if p\ <^ p2, the error of the pressure can be ar- 
bitrarily large. Note that the existence of this error 
does not imply the inconsistency of SPH. In this limit 
of ft, — > 0, the volume of the regime affected by this er- 
ror approaches to zero, which means the original differ- 
ential equation is restored almost everywhere. In other 
words, SPH satisfies the weak form of the original equa- 
tion. However, it means the convergence is slow and first 
order. 

One might think that this error is caused by an inade- 
quate initial thermal energy distribution. However, it is 
not the case. Even if we initialize the internal energy of 
particles near the contact discontinuity so that pressure 
is smooth, the particle distribution changes and discon- 
tinuity is regenerated. We, thus, need continuous adjust- 
ment to suppress the pressure error through out the time 
integration. Price's artificial conductivity (jPricd [2"008| ) 
provides such a continuous adjustment. Though the ar- 
tificial conductivity works beautifully in test calculations 
for the Kelvin-Helmholtz instability, whether its use in 
actual astrophysical simulation is justified or not is a bit 
questionable. First, in the case of the discontinuity of 
chemical composition, not only the jump in the internal 
energy but also that in the chemical composition should 
be smoothed but that is clearly not adequate. Second, 
the artificial heat conduction can significantly enhance 
the thermal relaxation of the system, which is again un- 
wanted. 

3. A DENSITY INDEPENDENT FORMULATION OF SPH 

In ^ we have seen that the standard SPH breaks down 
at the contact discontinuity because the continuity and 
differentiability of the density is necessary to guarantee 
the convergence of SPH approximation. The basic rea- 
son for this problem is the use of mj/pj as the volume 
element. Thus, if we use something else as the volume 
element, we might be able to avoid this difficulty alto- 
gether. 

3.1. Concept 

Here, we propose an alternative formulation of SPH in 
which we discretize Eq. [5] using the EOS of fluid, not the 
mass density. The new volume element is 



dr' 



(7 - l)mjUj 
Pi 



(24) 



Substituting Eq. [M] into Eq. [5] and using the gather 
summation, we obtain a new SPH approximation of 
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Fig. 1. — Density and pressure fields evaluated with the standard SPH and our new SPH around the contact discontinuity with the 
density ratio of 1 : 8. Equal mass particles are used for the first and second configurations (the left and middle columns). The positions 
of the less dense region is determined by taking the center of mass of the eight particles in the cube of the particle separation. The equal 
separation is used in the last configuration (the right column). In this configuration, the mass of particles in the less dense region is 1/8 of 
that of particles in the dense region. For the left and right columns, the constant neighbor number, 32 ± 2, is used. In the middle column, 
a constant kernel size of 0.03125 is used. The top row sho ws t he distribution of particles projected on the x — y plane. The second row 
shows the density of each SPH particle evaluated with Eq. 1161 The third row shows the density contrast between the evaluated density 
and true one. The fourth row shows corresponding pressure. In the bottom row, the pressure of each particle calculated with our new SPH 
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smoothed /: 

(/>(n): 



3 



(7-1) 



(25) 
(26) 



where Ui = rriiUi is the internal energy of particle i. By 
substituting / with the energy density, q = pu, we have 



3 



U,W{n,,h,). 



(27) 



where we used qi = {q){r). The gradient of (/) is given 

by 

(V/)(n)^^C/,^Viy(r,,,/i,)- (28) 

3 

We adopt Eqs. [53] and [57] as the basis of our new for- 
mulation. We derive the equations of motion and energy 
from this basis. We note that our new SPH is also Galilei 
invariant. 

One might think that the use of U for the calculation 
of the volume element would cause some inconsistency, 
since U is not a conserved quantity. The mass of a par- 
ticle is constant, and thus looks safer. In the following, 
we show that we can construct a consistent set of equa- 
tions using [/, and that it has many advantages over the 
standard SPH and that it retains important characteris- 
tics such as force symmetry and energy conservation. We 
first derive the energy equation and then equation of mo- 
tion. We then discuss the formulation for the estimate 
of the density and the implementation of the artificial 
viscosity. 

3.2. Energy Equation 

We need an expression of V • to derive the energy 
equation. We start with 

V{qv) = Vqv + qV ■ V, (29) 

which is obtained by replacing p in Eq. [T^] by q. By 
applying Eq. [28] to Eq [29] we obtain an analogy of Eq. 

q^V ■ Vi 



3 



The energy equation is then given by 



dui 



T.U3 



Pi_ 
Piqi 



(30) 



(31) 



Equation [31] contains pi since u is the energy per unit 
mass. The equation for Ui is obtained by multiplying 
Eq. [31] by mi: 

dUi m, \ - UiP^ 
IT 



Pi 



3 



-Vij ■ VWij. 



(32) 



Here, rui/pi is the volume associated with particle i 
which can be replaced by Ui/qi. Thus, we have 



dU, 



where we used P = (7 — \)q. 



(33) 



3.3. Equation of Motion 

From the energy equation, Eq. I33[ we derive the equa- 
tion of motion. The change in the internal energy of 
particles i and j due to their relative motion is 

^ + ^ = (7 - (- + -) v., ■ Vm,. (34) 



dt 



dt 



This change is the same as the change of the kinetic 
energy of particles with an opposite sign. Thus, we have 



/ dv^ 


dVj\ 


fdU, 




\ dt 


dt I 


\ dt ~ 


dt / 



mi + mj 

Substituting Eq. [34] into Eq. [35] we obtain 



dvi dvj 
~dt dt 



mi + mj f 1 1 

-7-1 '-U^Uj - + — 

m.imj \qi qj 



. (36) 

Since the motion of the center of mass of two particles is 
unchanged by the interaction of two particle, we have 



— (m^Dj + mjVj) = 0. 



(37) 



Thus, we have 



= -(7 - W^UJ (- + -) yW,j, (38) 
dt V* IjJ 

as the contribution of particle j to the equation of motion 
of particle i. 

The equation of motion for particle i is obtained by 
taking summation over neighbor particles: 



dvi 
' dt 



-(7-i)Ew,(i + l 



(39) 



The right-hand side of Eq. [39] contains only the energy 
U and energy density q. Thus, as far as q is smooth, Eq. 
[551 is likely to be well-behaved. The equation of motion 
of the standard SPH (Eq. [9]) requires both P and p 
are smooth. Thus, in our formulation, there is nothing 
special about the contact discontinuity. We can therefore 
expect that the treatment of the contact discontinuity is 
improved. We will see this in ij3.5l 

Note that Eq. [39] is mathemati cally equivalent to 
the eq uation of motion obtained by [Ritchie &: Thomas! 
(|2001|). while the deviation is completely different. 
[Ritchie fc Thomas! ([2001I ) started from Eq. [27] and den- 
sity estimate p = mq/U, but still tried to use standard 
SPH estimate of Eq. [B] In order to eliminate p from 
equation of motion, they used the following formal rela- 
tionship 

VP VP P„^ 

= + -V1, (40) 

P P P 



and formal identity 

VI = ^mj— VM^(rij,/i) ~ 0. 
Pj 



(41) 



Thus, their deviation was a heuristic modification of the 
standard SPH and they did not employ the volume el- 
ement U/q explicitly. We have shown that by choosing 
U/q as the volume element, we can derive a consistent 
set of SPH equations naturally. 
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3.4. Artificial Viscosity 

To deal with shocks, the standard SPH needs an ar- 
tificial viscosity term. Our new formulation also needs 
an artificial viscosity term. We utilize artificial viscos- 
ity terms which are widely used in simulations with the 
standard SPH. 

The viscosity term for the equation of motion is 

tfr. 

mj-^ = -m, ^ mj-Hij VWij, (42) 

3 

and the corresponding form of it for the energy equation 
is 



dt 



(43) 



where Ily is the function of the strength of the artificial 
viscosity. 

There are two types of artificial viscosity term, Ilij, 
whic h are commonly used . The most commonly used 
one (|Lattanzio et al.|[l985[ ) is 



n,; 



Pii 







Vi^ ■ r.ij > 0, 



(44) 



where a and /3 are the control parameters for the strength 
of the artificial viscosity, dj is the arithmetic average of 
the sound speeds of particles i and j, pij = 0.5{pi + Pj), 
and 

(45) 



r2. 



1^ 



The constant e is introduced to avoid the divergence and 
its fiducial value is ~ 0.01. 

The other one is proposed by iMonaghanl ()1997t ) from 
the analogy of the Riemann solver: 



n, 







(46) 



where 



Ci + Cj — 3wij and Wij = Vij ■ rij/rij. 



Since we have the density estimate p = q/u, we have 

(47) 



However, this modification of pij leads to unstable behav- 
ior under strong shocks. It seems to be safe to use the 
smoothed densities of particles i and j evaluated using 
Eq. [71 A consistent derivation of the artificial viscosity 
term in our scheme will be investigated in a forthcoming 

paper. 

We use the standard Balsara switch ()Balsaralll995f ) to 
suppress the shear viscosity. It is given by 



IV • V,; 



V -Vil + lV X Vi\ +ebCi/hi 



(48) 



and Hij^Baisara = 0-5{Fi + Fj)Ilij . Here eb is a small value 
(typically lO"*^). The rotation of velocity in the sta ndard 
SPH is found in literature (e.g.. lMonaghanlll992l ). The 
rotation of velocity in our SPH is calculated as follows: 



V X t)i ~ — ^ UjVtj X WW{r,j,hi). 



(49) 



3.5. Pressure in Contact Discontinuities 

The pressure around the contact discontinuity calcu- 
lated with our SPH equation is shown in the bottom 
panels of figure [T] In the case of the equal-mass parti- 
cle and the fixed neighbor number (the left panel), we 
can see that the jump of the pressure at the contact dis- 
continuity in our SPH is much smaller than that in the 
standard SPH. In the case of the constant kernel size (the 
middle panel), the result of our SPH is almost flat, while 
that of the standard SPH has a large error. 

In these two equal-mass cases, pressure still has small 
jumps at the contact discontinuity. The reason is that in 
both cases the distribution of particles is asymmetric. In 
the high-density region, the particle separation is smaller, 
resulting in small integration error. As a result, small 
error appears when the kernel contains the contribution 
from both low- and high-density regions. In the case 
of the equal separation of particles, there is no jump in 
the pressure distribution at the contact discontinuity, as 
shown in the rightmost panel. 

4. NUMERICAL EXPERIMENTS 

In this section, we show the results of several standard 
tests for fluid schemes, for both the standard SPH and 
our new SPH. In H4.1[ we describe our numerical code 
briefly. In ^4.21 we show the results of the shock tube 
tests. Then we show the evolution of system which is 
initially in hydrostatic equilibrium in ^4.31 In ^AA\ and 
i j4.5i tests for two important fluid instabilities are carried 
out. Finally, we show the results of the bl ob tests which 
was first proposed bv lAgertz et al.l ()2007[ ). In all tests, 
our new SPH shows much better result compared to that 
of the standard SPH. 

4.1. Numerical Method 

We used ASURA, a parallel 7V-body/SPH code, as a 
framework of current numerical experiments. ASURA 
adopts the leap-frog method for the time- integration. For 
simplicity, we used the shared steps with variable time- 
steps. The time-step is given by 

dt = mindti, (50) 



where 



dt, = a 



2h,, 



CFL- 



(51) 



and CcFL = 0.3. 

For the standard SPH, we first evaluated the densi- 
ties of particles using Eq. [161 Then, we calculated the 
pressure gradient and the time-derivative of the internal 
energy using Eqs. [T7land[T51 In our new SPH, we eval- 
uated q using Eq. [27l first, and then we calculated the 
pressure gradient and the time derivative of the internal 
energy using Eqs. [39land[33l We used Eq. [46las the arti- 
ficial viscosity term in both cases and we adopted a = 1. 
The Balsara switch was also applied. To avoid the tensile 
instability, we us ed a first derivative of the ke rnel which 
has a cuspy core (|Th omas & Couchmanl[l99l ). 

The kernel size of each particle is determine to keep the 
number of neighbor particles within the range of 32 ± 2. 
As an exception, in the one dimensional tests shown in 
the kernel size is evaluated by 



P 



(52) 
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where rj = 1.2 for the Sod's shock tube tests and rj = 2.4 
for the strong shock tube tests. 

For brevity, we sometimes express the standard SPH 
as SSPH. 

4.2. Shock Tube Tests 

The Sod shock tube ()Sodlll97"8[ ) is the most basic test 
for numerical schemes for compressible fluid. This test 
shows the shock capturing ability of schemes. In SPH, 
not only the profile of the shock front but also the be- 
havior of the contact discontinuity is important. Here, 
we show the results of one- and three-dimensional shock 
tube tests. 

The setup is as follows. We prepared the periodic do- 
main of —1 < X < 1 for the one-dimensional tests and 
-1 < a; < 1, -1/16 <y < 1/16, and -1/16 < z < 1/16 
for the three dimensional tests. The initial condition is 
given by 

(p=l,P=l,v = 0, x<0, 

\p = 0.25, P = 0.1795, w = 0, a; > 0. ^ ' 

To express this initial condition, we use equal-mass parti- 
cles and place 800 and 200 particles in the left and right 
domains, respectively, regularly in the one-dimensional 
tests. In three dimensional tests, we place 40000 and 
10000 particles in the left and right domains, respec- 
tively, and a glass-like particle distribution was used. We 
set 7 = 1.4 and gave the internal energy to each particle 
to ensure the given P. 

In addition to the Sod shock tube, we performed a one- 
dimensional strong shock test. The initial condition for 
this test is given by 

fp = 1,P = 1000,u = 0, a; < 0, , , 

\p = 1,P = 0.01,w = 0, a;>0. ^ ' 

We use 1000 equal-mass particles in the computational 
domain of — 1 < a; < 1 with the equal separation. 

Figure [2] shows the results of the one-dimensional shock 
tube tests with the standard SPH and our new SPH. The 
density (upper row) and pressure (bottom row) of each 
particle are plotted by circles. The red curves represent 
the analytic solutions. 

The standard SPH reproduce the analytic solution of 
the density distribution well. The shock front is resolved 
by ~ 7 particles. The jump of the density at the contact 
discontinuity is resolved by a similar number of particles. 
The pressure shows large variations near the contact dis- 
continuity, though it should be constant. Since Eq. IHlof 
the standard SPH contains a large error near the con- 
tact discontinuity, in order to achieve zero acceleration, 
pressures of particles must have large variations. This 
result is the same as th e results of previous work s with 
the standard SPH fe.g.. lST)ringell[20n5l : IPricai2n08l\ 

In our new SPH, unlike the case of the standard SPH, 
the pressure around the contact discontinuity does not 
show a large jump. The reason is simply that the energy 
density is used instead of the mass density. The energy 
density is constant at the contact discontinuity. The rea- 
son why there is a small change in the pressure is that 



T 




Q I \ \ \ \ \ I I \ \ \ \ \ 

-0.4-0.2 0.2 0.4 -0.4-0.2 0.2 0.4 

Position Position 



T 




Q I \ \ \ \ \ I I \ \ \ \ \ 

-0.4-0.2 0.2 0.4 -0.4-0.2 0.2 0.4 



Position Position 

Fig. 2. — The results of the one dimensional shock tube tests 
for SSPH and our new SPH at t = 0.1. Density (upper row) and 
pressure (bottom row) are shown. Circles indicate the physical 
quantities of each SPH particle, while red curves represent the 
analytic solutions. Insets in the pressure panels arc the close-up 
views around the contact discontinuity. 

the particle separation changes at the contact disconti- 
nuity. As we showed in figure [U our new SPH still has 
small error in the pressure, due to the finite number of 
particles in the kernel. This error caused the change in 
the pressure in figure [21 

The results of the three dimensional shock tube tests 
for the standard SPH and our SPH are shown in figure 
[21 In this figure, the circles represent average values of 
particles in bins with the width of the mean particle sep- 
aration at the high density part. Again, we can see a 
variation in the pressure around the contact discontinu- 
ity in the case of the standard SPH. In the case of our 
SPH, there is no such variation. 

Figure [31 shows the results of the strong shock tube 
tests for the standard SPH and our new SPH. The shock 
front and the contact discontinuity in the density dis- 
tribution is well reproduced in the both cases. In this 
extreme test, both runs show jumps in the pressure dis- 
tribution around the contact discontinuity. The absolute 
value of the pressure jump in our SPH is much smaller 
than that in the standard SPH. The jump found in the 
pressure in our SPH is caused by the asymmetry in the 
particle distribution (see Overall, our SPH can 

handle such a strong shock problem, even when a very 
large pressure jump exists initially. This result is quite 
reassuring. In our new SPH, it is assumed that pressure 
is continuous, which is not a valid assumption at the 
shock front. Thus, it could fail to capture very strong 
shocks. The result shown in figure [31 shows that is not 
the case and our new SPH can handle very strong shocks. 

4.3. Hydrostatic Equilibrium Tests 

As is shown in fj2l in the standard SPH, particles feel 
unphysical repulsive force at the interface of the contact 
discontinuity. Therefore, in order to establish the hydro- 
static equilibrium, the distance between particles at the 
different sides of the contact discontinuity must become 
larger than the "true" value. What is the consequence of 
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Fig. 3. — The results of the three dimensional shock tube tests 
for SSPH and our new SPH at t = 0.1. Density (upper row) and 
pressure (bottom row) are shown. Circles indicate the averaged 
physical quantities of SPH particles, while red curves represent the 
analytic solutions. 
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Fig. 4. — The same as figure (2] but for the strong shock tube 
tests at t = 0.012. 



this repulsive force? Here, we show the result of a sim- 
ple test which helps us to understand the problem of the 
unphysical repu l sive fo rce. Similar test has been used in 
iHeB fc Sprineeil (IMol) . 

We follow the evolution of two fluids with different val- 
ues of density, but with the same pressure. We performed 
two-dimensional tests. The computational domain is a 
square of the unit size, < a: < 1 and < y < 1, with a 
periodic boundary condition. Initial conditions are 









{: 


p = 


2.5, 


7 = 


5/3 



0.25 < a; < 0.75 and 0.25 <y < 0.75, 
otherwise, 



(56) 

(57) 
(58) 

We tried two different realizations. In the first one, the 



particle mass is the same for the entire computational 
region. Thus, the inter-particle distance is smaller in the 
high density region. In the second one, particles in the 
high density region is four times more massive than par- 
ticles in the low-density region. In both cases, particles 
arc initially in a regular grid. For the equal-mass case, 
the number of particles in the dense region is 4096 and 
that in the ambient is 3072. For the equal-separation 
cases, those arc 1024 and 3072, respectively. Initial ve- 
locities of particles were set to zero. Since the system is 
initially in the hydrostatic equilibrium, particles should 
not move, except for small local adjustments. 

Figure [5] shows the time evolution up to t = 8. There 
is a clear difference between the result of the standard 
SPH and that of our SPH. With the standard SPH, the 
high-density region, which initially has a square shape, 
quickly becomes rounder and almost completely circular 
by i = 8. We can understand this unphysical round- 
ing as follows. As we stated in fj2] and §3.51 unphysical 
repulsive force between particles operates at the contact 
discontinuity. We can see the effect of this force in the de- 
velopment of the gap of the distribution of particles near 
the boundary of two fluids. Because of this gap, the bulk 
of the system is slightly compressed. The system seeks to 
achieve the energy minimum, by minimizing the surface 
area of the contact discontinuity. Thus, the high-density 
region evolves to a circular shape, which minimizes the 
length of the boundary. In other words, the repulsive 
force effectively adds the "surface tension" . 

Our new SPH gives a far better solution, as we can 
see in the lower two rows of figure [S] The overall square 
shape remains there till the end of the simulation in the 
equal-mass case. The result of the unequal-mass case is 
even better. The equation of motion of our SPH elimi- 
nates the unphysical surface tension completely. 

Figure [6] shows the final state of the two-fluid system 
with the density contrast of 64. Our SPH handles the 
system without any problem (right panel). On the other 
hand, in the calculation with the standard SPH, a wide 
and empty ring structure is formed between two fluids. 

4.4. Kelvin- Helmholtz Instability Tests 

After the work bv lAgertz et al.l ()2007() which demon- 
strated clearly that the standard SPH cannot deal with 
the Kelvin-Helmholtz instability correctly, many re- 
searchers proposed modifications of SPH to solve the 
problem (see fjl])- In this section, we investigate how 
our new SPH handles the Kelvin-Helmholtz instability. 

We prepared a two-dimensional computational do- 
main, < a; < 1 and < y < 1. The periodic boundary 
condition was used. The density is 



l{=pi) 0<y < 0.25, 0.75 < 2/ < 1, 
2(= ph) 0.25 <y < 0.75. 



(59) 



We used equal-mass particles. The numbers of particles 
in the high and low dense regions are 131072 and 65522, 
respectively. We set P = 2.5 and 7 = 5/3. The high and 
low density regions had the initial velocities of Vx,h ~ 0.5 
and Vx,i = —0.5 in the x direction, respectively. 

We have used A^nb = 32 ± 2 as the neighbor num- 
ber. This value might seem a bit large, but we found 
it guarantees the good sampling of the particles in the 
low-density region at the interface. When we used 
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Fig. 5. — Snapshots of a two-fluid system at t = 0.1, 0.3, 0.5, 1 and 8. The red and blue points indicate the positions of particles with 
p = 4 and p = 1, respectively. The upper two rows are the results of SSPH. while the lower two rows are those of our SPH. The first and 
third rows show the results of the equal-mass cases, whereas the second and fourth rows show those of the equal separation and unequal 
mass cases. 



■SSPH & Mass ratio 1 :64:Sv;;S 
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Fig. 6. — The final state (t = 8) of a two fiuid system with the 
density contrast of 64. The red and blue points are the positions of 
particles with p = 64 and 1, respectively. The particle separation 
is constant and the particle mass difference is 1:64. 



A'lib = 16 ± 2, the variation of the pressure at the in- 
terface becomes too large. For the artificial viscosity, we 
used a = 1 with the Balsara switch. 

We added a small velocity p erturb ation to the particles 
near the interfaces, following iPricd ([20081 ). The velocity 
perturbation in the y direction is as follows: 



Aw,, 



Asm[~2Tr{x + 0.5)/X], \y - 
Asin[27r(a--|-0.5)/A], \y - 



where A = 1/6 and A = 0.025. 



0.25| < 0.025 
0.75| < 0.025, 
(60) 



The time-scale of the growth of the Kelvin-Helmholtz 
instability is 



Tkh 



MPh+Pi) 



(61) 



For our test setup, r^h = 0.35. We followed the evolution 
up to t = Srich- 

The results are shown in figure [71 The difference be- 
tween two results is clear. In the run with the stan- 
dard SPH, perturbations grow till t = Tkh, but the un- 
physical surface tension inhibited the growth of roll-like 
structures. The stretched high-density fluids break apart 
{t ~ 4Tkh) and form blobs {t = Sr^h)- This evolu- 
tion is com pletely different from those obtained by Euler 
codes fe.g.. lAgertz et aLll2007l) . On the other hand, our 
new SPH shows a very good result which is compara- 
ble to those with Euler codes and with SPH with the 
iRitchie fc Thomai ()2001l) equation o f motion or the ar- 
tificial conductivity (see lPricel[2008[) . 

Figure [5] shows the cross section of the pressure distri- 
bution along the j/-axis. We can see that a very large 
pressure jump exists around the contact interfaces, in 
the case of the standard SPH. The surface tension at the 
interface of the two fluids prevents the normal growth 
of the Kelvin-Helmholtz instability. On the other hand, 
there is no such jump in the case of our SPH. Since the 
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Fig. 7. — The density maps from the two dimensional shear flow test at t = 1, 2, 4 and 8 T\h- The upper panels show the results of SSPH, 
while the bottom panels show those of our SPH. The color code of the density is given at the bottom. 



pressure and particle distribution is well-behaved at the 
interface, the growth of the Kelvin-Helmholtz instability 
is not suppressed. 

4.5. Ray leigh- Taylor Instability Tests 

lAbell (|2011[ ) demonstrated that the standard SPH can- 
not follow the development of the Ray leigh- Taylor insta- 
bility correctly. We show the result with our SPH as well 
as that with the standard SPH. 

The initial setup is as follows. We prepared the two 
dimensional computational domain of < x < 1 and 
< y < 1. We placed two fluids separated at y = 0.5. 
The density just above (below) the interface was set to 
Ph = 2 (pi = 1). These two fluids were initially in the 
hydrostatic equilibrium. Further, we assumed that each 
fluid was initially isoentropic. The density distributions 
of these fluids in the vertical direction arc given by 



Pi 



Ph 



1 



7-1 pig(y-0-5) 

7 Po 



1 I 7-1 Phg{y-0-5) 

Po 



y < 0.5, 
y > 0.5, 



(62) 



where g = —0.5 is the gravitational constant, Pq = 10/7 
is the value of pressure at the interface, and 7 = 1.4. The 
initial density and entropy profiles arc shown in figure [9l 
To ensure the initial density distribution given by Eq. 1621 
we first placed equal-mass particles on the regular grid 
with the separation of 1/512. Then, we adjusted the 
vertical separation of each particle set having the same y 
to reproduce the density distribution. The particle mass 
was set to 5.7x 10~^ and the total number of particles was 
247296. The periodic boundary condition was imposed 
on the X direction. Particles with y < 0.1 and y > 0.9 
were fixed at the initial positions and they were not allow 
to change their internal energy. 

The velocity perturbation in the vertical direction was 
added as the seed of the instabilities. We carried out 



runs with two kinds of the seed. For the first test, we 
added the velocity perturbation to particles in the range 
of 0.3 < y < 0.7, and the form of the perturbation is 

Avy{x,y) ^ Syy[l + cos{4:TTx)]{l + cos[5TT{y~0.5)]}. (63) 

We set Svy = 0.025. For the second test, we added the 
velocity perturbation of the form: 

Avy{x,y) 



40 
j=20 



cos{kjx) exp{—Q.05kj\y 

kj 



and 



Ph - Pi 
Ph + Pi 



0.5|), 
(64) 

(65) 



where rij is the linear growth rate of the Ray leigh- Taylor 
instability, and kj = ji: / L{= 1) is the wave number of 
the perturbation. The amplitude of each mode, Oj, was 
drawn from a Gaussian distribution with the variance 
of unity at random. Th is initial velocity perturbation is 
based on lYoung^ ()1984[ ) with slight modifications. Veloc- 
ities of the particles outside the perturbed region was set 
to zero. We call these two tests single-mode and multi- 
mode tests, respectively. 

In figure [TUl the growth of the Rayleigh- Taylor insta- 
bility in the case of the single-mode test is shown. The 
Rayleigh- Taylor instability develops in calculations with 
both of the standard SPH and our SPH, but the struc- 
tures of them arc quite different. The unphysical surface 
tension of the standard SPH again prevents the develop- 
ment of the fine structures on the surface of the Rayleigh- 
Taylor fingers. Thus, the result looks quite different from 
those obtained with Euler schemes. On the other hand, 
in the calculation with our SPH, the overall evolution of 
the Rayleigh- Taylor instability in our SPH shows excel- 
lent agreement with thos e with Euler sc hemes and the 
moving mesh scheme (see lSpringelll2010al) . 
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Fig. 8. — Pressure of each particle along the y direction at t = 0.4 t^^j. The left panel shows the result of SSPH, whereas the right panel 
shows that of our SPH. Particles initially in the high- (low-) density region are expressed with red (blue) points. 
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Fig. 9. — Initial distributions of density and entropy in the 
vertical direction. Solid and dotted curves indicate density and 
entropy, respectively. 



Figure [TT] shows the growth of the Raylcigh- Taylor 
instability with the multi-mode perturbations with the 
standard SPH and our new SPH. The global phase mix- 
ing of fluids can be seen in the result with our SPH. On 
the other hand, due to the unphysical surface tension, the 
mixing is significantly suppressed in that of the standard 
SPH. The distribution of two fluids looks like a mixture 
of oil and water. 

4.6. Blob Tests 

As the final test, we performed the blob test proposed 
bv lAgertz et al.l (|2007( ). This test incorporates both the 
Kelvin-Helmholtz and Rayleigh- Taylor instabilities. 

We used Read ' s initial condition of the blob test 
(iRead et al.l[2010l: iRead fc Havfieldl IMl 2. The com- 
putational domain was < x < 2000 kpc, < y < 

^ We obtained the initial condition from the following URL: 
|\proteclj^ittp://www-theorie.physik.uzh.ch/astrosim/code/| 



2000 kpc, and < z < 6000 kpc, and the periodic bound- 
ary condition was imposed. A cold cloud of the density 
= 3.13 X 10~^ in the mass unit of 2.3 x lO'^Afoand 
the length unit of 1 kpc and temperature Tc = lO^K was 
centered at {x,y,z) = (1000 kpc, 1000 kpc, 2000 kpc). 
The radius of this cloud was 197 kpc. This cloud was 
embedded in the diffuse ambient gas of which density 
and temperature were pa = 3.13 x 10~* and Ta = 10^ K, 
respectively. The ambient gas had the velocity of Vz = 
1000 kms^^. Thus, the Mach number of the flow to the 
cloud was 2.7. The total number of particle for the sys- 
tem is 4643283. We integrated the system up to t = 5rkh, 
where r^h = 2 Gyr is the typical growth tim e-scale of the 
Kelvi n- Hclmholtz instability in this test (jAgertz et al.l 
[2OO7I ). 

Figure [12] shows the snapshots of the cloud core. The 
upper and lower panels are the results with the stan- 
dard SPH and our SPH, respectively. Their evolutions 
were quite different. The blob simulated with the stan- 
dard SPH retained the single cloud structure until the 
late stage of the simulation. This behavior is consistent 
with t hose with the standard SPH shown in lAgertz et al.l 
(|2007( ). In contrast, the blob surface was disrupted in the 
run with our SPH, due to the growth of the instabilities 
on the surface. The blob fragmented into several peaces 
and mixed eventually with the ambient gas. This behav- 
ior is similar to those obtained by Euler codes. 

The evolution of the blob mass is shown in figure [T3l 
Here we show t he mass of gas with p > 0.64 and T < 
0.9 Ta, following lAgertz et al.l (|2007t ). At t = 2.5 Tkh, the 
blob mass in the run of our SPH became ^ 10 % of the 
initial mass. This result is consist ent with the result s 
of the Euler codes (see figure 6 in lAgertz et al.l l2007t ) . 
The evolution of the blob mass in the standard SPH was 
much slower compared to that in our SPH. 

5. SUMMARY AND DISCUSSION 

In this paper, we described an alternative formulation 
of SPH in which the energy density is used as the basis of 
the smoothing instead of the mass density. In our formu- 
lation, the mass of particles is not used in the evaluation 
of the right-hand sides of the energy equation and the 
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Fig. 10. — The density maps of the two dimensional Raylcigh- Taylor instability tests at t = 0.5, 3, 4 and 5. The upper panels show the 
results of SSPH, while the bottom panels show those of our SPH. The color code of the density is given at the bottom. 



equation of motion. As a result, the large error of force 
estimate at the contact discontinuity, which is unavoid- 
able with the standard SPH, disappears completely in 
our SPH. Not surprisingly, our SPH can handle contact 
discontinuities and the Kelvin-Helmholtz and Rayleigh- 
Taylor instabilities without difficulty. The behavior of 
the shock in our new SPH is essentially the same as 
that in the standard SPH. Since the equations used in 
our SPH are almost identical to those of in the stan- 
dard SPH except that energy density is used in place of 
mass density p. Modification of existing SPH code to use 
our scheme is simple and straightforward. In particular, 
there is no increase in the calculation cost. Equations 
which are not derived in this paper, such as the diffusion 
e quatio n ([Bro okshaw 1985), can be derived easily. 

iPricd ()2008D improved the treatment of the Kclvin- 
Helmholtz instability of the standard SPH, by applying 
artificial conductivity at the contact discontinuity. Un- 
like the artificial viscosity, artificial conductivity intro- 
duces the dissipation not in the original set of equation. 
Our SPH docs not need such additional dissipation, and 
thus the contact discontinuity is kept sharp. 



One might think our result contradicts with the re- 
quirement that a ll quantities in SPH must be smooth 
(|Monaghanl Il997t ). However, it is obvious that in our 
SPH; all quantities in the right-hand side of the equa- 
tions are smooth, since the only discontinuous quantity 
is the density and it does not appear in the right-hand 
sides. Thus, our results does not contradict with Mon- 
aghan's requirement. 

In this paper, we discuss the treatment of ideal gas only. 
We arc currently working on the extension to non-ideal 
fluid, and the result will be given in the forthcoming 
paper. 



Some of the numerical tests were carried out on the 
Cray XT4 system in the Center for Computational As- 
trophysics at the National Astronomical Observatory of 
Japan. This work is supported by HPCI Strategic Pro- 
gram Field 5 'The origin of matter and the universe' and 
Grant-in- Aid for Scientific Research (21244020) of Japan 
Society for the Promotion of Science, Ministry of Educa- 
tion, Culture, Sports, Science and Technology, Japan. 
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Fig. 12. — The density maps at t = 0.25, 1.0, 1.75 and 2.5 r]jii. The upper and lower panels show the results with SSPH and our SPH, 
respectively. The color code of the density is given at the bottom. 




